We will analyse the data collected by Jones (Unpublished BSc dissertation, University of Southampton, 1975). The aim of the study was to define if the probability of having Bronchitis is influenced by smoking and/or pollution.
The data are stored under data/Bronchitis.csv and contains information on 212 participants.
Lets starts by
read.csv()Bronchitis = read.csv("data/Bronchitis.csv",header=TRUE)
plot(Bronchitis$cigs,Bronchitis$bron,col="blue4",
ylab = "Absence/Presense of Bronchitis", xlab = "Daily number of cigarettes")
abline(h=c(0,1),col="light blue")Lets
glm() and by means of the function gamlss() of the library gamlss.glm function : Use the function summary() to display the results of an R object of class glm.fit.glm = glm(bron~cigs,data=Bronchitis,family=binomial)
library(gamlss)
fit.gamlss = gamlss(bron~cigs,data=Bronchitis,family=BI)## GAMLSS-RS iteration 1: Global Deviance = 181.7072
## GAMLSS-RS iteration 2: Global Deviance = 181.7072
##
## Call:
## glm(formula = bron ~ cigs, family = binomial, data = Bronchitis)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4418 -0.5472 -0.4653 -0.4405 2.1822
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.2840 0.2731 -8.365 < 2e-16 ***
## cigs 0.2094 0.0376 5.567 2.59e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 221.78 on 211 degrees of freedom
## Residual deviance: 181.71 on 210 degrees of freedom
## AIC: 185.71
##
## Number of Fisher Scoring iterations: 4
Let’s now define the estimated probability of having bronchitis for any number of daily smoked cigarette and display the corresponding logistic curve on a plot:
plot(Bronchitis$cigs,Bronchitis$bron,col="blue4",
ylab = "Absence/Presense of Bronchitis", xlab = "Daily number of cigarettes")
abline(h=c(0,1),col="light blue")
axe.x = seq(0,40,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])/(1+exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2]))
lines(axe.x,f.x,col="pink2",lwd=2)As for linear models, model selection may be done by means of the function anova() used on the glm object of interest.
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: bron
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 211 221.78
## cigs 1 40.07 210 181.71 2.45e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Lets assess is the model fit seems satisfactory by means
plot() on an object of class glm,plot() on an object of class gamlss,## GAMLSS-RS iteration 1: Global Deviance = 181.7072
## GAMLSS-RS iteration 2: Global Deviance = 181.7072
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = -0.05227133
## variance = 0.9386897
## coef. of skewness = 0.02624943
## coef. of kurtosis = 2.950084
## Filliben correlation coefficient = 0.9984055
## ******************************************************************
# long format:
long = data.frame(mi = rep(c("MI","No MI"),c(104+189,11037+11034)),
treatment = rep(c("Aspirin","Placebo","Aspirin","Placebo"),c(104,189,11037,11034)))
# short format: 2 by 2 table
table2by2 = table(long$treatment,long$mi)
#
chisq.test(table2by2)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table2by2
## X-squared = 23.782, df = 1, p-value = 1.079e-06
##
## 2-sample test for equality of proportions with continuity correction
##
## data: table2by2[, "MI"] out of apply(table2by2, 1, sum)
## X-squared = 23.782, df = 1, p-value = 1.079e-06
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## -0.010570828 -0.004440228
## sample estimates:
## prop 1 prop 2
## 0.009334889 0.016840417
##
## Call:
## glm(formula = I(mi == "MI") ~ treatment, family = "binomial",
## data = long)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.1843 -0.1843 -0.1370 -0.1370 3.0574
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.66462 0.09852 -47.348 < 2e-16 ***
## treatmentPlacebo 0.59763 0.12283 4.865 1.14e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 3122.5 on 22363 degrees of freedom
## Residual deviance: 3097.8 on 22362 degrees of freedom
## AIC: 3101.8
##
## Number of Fisher Scoring iterations: 7
The dataset students.csv shows the number of high school students diagnosed with an infectious disease for each day from the initial disease outbreak.
Lets
read.csv()cases) as a function of the days since the outbreak (variable day).students = read.csv("data/students.csv",header=TRUE)
plot(students$day,students$cases,col="blue4",
ylab = "Number of diagnosed students", xlab = "Days since initial outbreak")
abline(h=c(0),col="light blue")Lets
glm() and by means of the function gamlss() of the library gamlss.glm function : Use the function summary() to display the results of an R object of class glm.fit.glm = glm(cases~day,data=students,family=poisson)
library(gamlss)
fit.gamlss = gamlss(cases~day,data=students,,family=PO)## GAMLSS-RS iteration 1: Global Deviance = 389.1082
## GAMLSS-RS iteration 2: Global Deviance = 389.1082
##
## Call:
## glm(formula = cases ~ day, family = poisson, data = students)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.00482 -0.85719 -0.09331 0.63969 1.73696
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.990235 0.083935 23.71 <2e-16 ***
## day -0.017463 0.001727 -10.11 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 215.36 on 108 degrees of freedom
## Residual deviance: 101.17 on 107 degrees of freedom
## AIC: 393.11
##
## Number of Fisher Scoring iterations: 5
Let’s now define the estimated probability of having bronchitis for any number of daily smoked cigarette and display the corresponding logistic curve on a plot:
plot(students$day,students$cases,col="blue4",
ylab = "Number of diagnosed students", xlab = "Days since initial outbreak")
abline(h=c(0),col="light blue")
axe.x = seq(0,120,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])
lines(axe.x,f.x,col="pink2",lwd=2)As for linear models, model selection may be done by means of the function anova() used on the glm object of interest.
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: cases
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 108 215.36
## day 1 114.18 107 101.17 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Lets assess is the model fit seems satisfactory by means
plot() on an object of class glm,plot() on an object of class gamlss,## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = 0.01029368
## variance = 0.895386
## coef. of skewness = -0.444238
## coef. of kurtosis = 2.638977
## Filliben correlation coefficient = 0.9881143
## ******************************************************************
Analyse further the Bronchitis data of Jones (1975) by
poll),cigs and poll.Lets plot the data first.
Bronchitis = read.csv("data/Bronchitis.csv",header=TRUE)
# plot
plot(Bronchitis$poll,Bronchitis$bron,col="blue4",
ylab = "Absence/Presense of Bronchitis", xlab = "Pollution level")
abline(h=c(0,1),col="light blue") No obvious relationship between pollution and bronchitis is visible by means of this plot.
Lets fit a model assuming that the probability of getting bronchitis is a function of the pollution level.
##
## Call:
## glm(formula = bron ~ poll, family = binomial, data = Bronchitis)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1018 -0.7254 -0.6004 -0.4944 2.0796
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.65578 2.59070 -3.341 0.000834 ***
## poll 0.12482 0.04341 2.876 0.004032 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 221.78 on 211 degrees of freedom
## Residual deviance: 213.26 on 210 degrees of freedom
## AIC: 217.26
##
## Number of Fisher Scoring iterations: 4
The intercept of the previous fit allows to define the probability of getting bronchitis when the level of pollution equals 0. As a zero level of pollution is (i) out of range (ii) not a realistic value, we will
poll_centered defined as the pollution level minus the mean so that the intercept corresponds to the probability of getting bronchitis for an average pollution level in Cardiff,# fit2:
Bronchitis$poll_centered = Bronchitis$poll-mean(Bronchitis$poll)
fit.glm = glm(bron~poll_centered,data=Bronchitis,family=binomial)
library(gamlss)
fit.gamlss = gamlss(bron~poll_centered,data=Bronchitis,family=BI)## GAMLSS-RS iteration 1: Global Deviance = 213.2586
## GAMLSS-RS iteration 2: Global Deviance = 213.2586
##
## Call:
## glm(formula = bron ~ poll_centered, family = binomial, data = Bronchitis)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.1018 -0.7254 -0.6004 -0.4944 2.0796
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.34808 0.17576 -7.670 1.72e-14 ***
## poll_centered 0.12482 0.04341 2.876 0.00403 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 221.78 on 211 degrees of freedom
## Residual deviance: 213.26 on 210 degrees of freedom
## AIC: 217.26
##
## Number of Fisher Scoring iterations: 4
Lets
## GAMLSS-RS iteration 1: Global Deviance = 213.2586
## GAMLSS-RS iteration 2: Global Deviance = 213.2586
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = 0.02469888
## variance = 0.8838803
## coef. of skewness = -0.1231663
## coef. of kurtosis = 2.542364
## Filliben correlation coefficient = 0.9962534
## ******************************************************************
## GAMLSS-RS iteration 1: Global Deviance = 213.2586
## GAMLSS-RS iteration 2: Global Deviance = 213.2586
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = -0.01887182
## variance = 1.227365
## coef. of skewness = -0.1119096
## coef. of kurtosis = 3.24346
## Filliben correlation coefficient = 0.9963054
## ******************************************************************
## GAMLSS-RS iteration 1: Global Deviance = 213.2586
## GAMLSS-RS iteration 2: Global Deviance = 213.2586
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = -0.07470508
## variance = 1.137648
## coef. of skewness = -0.04981677
## coef. of kurtosis = 2.853494
## Filliben correlation coefficient = 0.9982264
## ******************************************************************
# plot fit
plot(Bronchitis$poll_centered,Bronchitis$bron,col="blue4",
ylab = "Absence/Presense of Bronchitis", xlab = "Pollution level")
abline(h=c(0,1),col="light blue")
axe.x = seq(-10,10,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])/(1+exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2]))
lines(axe.x,f.x,col="pink2",lwd=2) Model check suggests a good fit. Lets finally check if the interaction is significant:
# interaction ?
fit.glm = glm(bron~poll_centered*cigs,data=Bronchitis,family=binomial)
summary(fit.glm)##
## Call:
## glm(formula = bron ~ poll_centered * cigs, family = binomial,
## data = Bronchitis)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4389 -0.5698 -0.4144 -0.2927 2.3970
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.394978 0.294856 -8.123 4.57e-16 ***
## poll_centered 0.158671 0.071836 2.209 0.0272 *
## cigs 0.214105 0.038699 5.533 3.16e-08 ***
## poll_centered:cigs -0.005292 0.010257 -0.516 0.6059
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 221.78 on 211 degrees of freedom
## Residual deviance: 173.95 on 208 degrees of freedom
## AIC: 181.95
##
## Number of Fisher Scoring iterations: 5
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: bron
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 211 221.78
## poll_centered 1 8.519 210 213.26 0.003515 **
## cigs 1 39.044 209 174.21 4.143e-10 ***
## poll_centered:cigs 1 0.264 208 173.95 0.607347
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Interaction is not significant.
The file myocardialinfarction.csv indicates if a participant had a myocardial infarction attack (variable infarction) as well the participant’s treatment (variable treatment).
Does Aspirin decrease the probability to have a myocardial infarction attack ?
Lets (i) import the dataset, (ii) change the levels of the factor treatment so that Placebo corresponds to the reference group, (iii) and finally plot the (sample) probabilities to get an attack by treatment group
# import
myocardialinfarction = read.csv("data/myocardialinfarction.csv")
# by default, Aspirin is the reference group as the alphabetic order is used
myocardialinfarction$treatment = factor(myocardialinfarction$treatment,
levels=c("Placebo","Aspirin"))
# plot
par(mfrow=c(1,1),mar=c(3,4,3,1))
pi.group = tapply(myocardialinfarction$infarction=="attack",myocardialinfarction$treatment,mean)
table.group = tapply(myocardialinfarction$infarction=="attack",myocardialinfarction$treatment,table)
temp = barplot(pi.group,plot=FALSE)
barplot(pi.group,ylab="Probability",xlab="",
main = "Probability of myocardial infarction\n per treatment group",names=rep("",2),
cex.axis=.6,axes=FALSE,cex.main=1.4)
axis(2,las=2,cex.axis=.8)
axis(1,temp[,1],names(pi.group),cex.axis=1.25,tick=FALSE)
for(gw in 1:length(pi.group)){
text(temp[gw],pi.group[gw]/2,.p(table.group[[gw]]["TRUE"]," / ",sum(table.group[[gw]])),
col="red",cex=1.5)
} The barplot seems to suggest that the treatment (aspirin) reduces the risk of myocardial infarction. Lets fit a logistic model to assess if this difference is significant. Note that, in this case (a dichotomous outcome and a dichotomous predictor), a test of equality of proportions or an independence test could also do the job. With a logistic model, other predictors could easily be added to the model and the beta parameter corresponding to the treatment can be interpreted by means of odd ratios (or relative risk ratios when prevalences are small, as we will note at the end of this practical).
# model fit
fit.glm = glm(I(infarction=="attack")~treatment,data=myocardialinfarction,family=binomial)
summary(fit.glm)##
## Call:
## glm(formula = I(infarction == "attack") ~ treatment, family = binomial,
## data = myocardialinfarction)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.1859 -0.1859 -0.1376 -0.1376 3.0544
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.04971 0.07337 -55.195 < 2e-16 ***
## treatmentAspirin -0.60544 0.12284 -4.929 8.28e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 3114.7 on 22070 degrees of freedom
## Residual deviance: 3089.3 on 22069 degrees of freedom
## AIC: 3093.3
##
## Number of Fisher Scoring iterations: 7
# test of equality of (independent) proportions
prop.test(unlist(lapply(table.group,function(x)x[2])),unlist(lapply(table.group,sum)))##
## 2-sample test for equality of proportions with continuity correction
##
## data: unlist(lapply(table.group, function(x) x[2])) out of unlist(lapply(table.group, sum))
## X-squared = 24.429, df = 1, p-value = 7.71e-07
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.004597134 0.010814914
## sample estimates:
## prop 1 prop 2
## 0.01712887 0.00942285
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: matrix(unlist(table.group), ncol = 2)
## X-squared = 24.429, df = 1, p-value = 7.71e-07
The three methods lead to the same conclusion: there is a significant difference between the probabilities of having a myocardial infarction of the two treatment groups. Note that the two last methods get exactly the same results (as they use the same test X-squared statsitic).
Lets define the fitted probabilities to get an attack and compare them to the sample probabilities: they should match:
## [1] 0.01712891
## Placebo
## 0.01712887
## [1] 0.009422852
## Aspirin
## 0.00942285
Finally, note that when prevalence are small, the exponential of the logistic regression corresponding to the treatment may also be interpreted as relative risk ratios. Indeed:
## [1] 0.5501137
## [1] 0.5458342
## [1] 0.4290414 0.6944202
Thus, aspirin strongly reduces the risk of myocardial infarction.
This data set is derived from Agresti (2007, Table 3.2, pp.76-77). It gives 6 variables for each of 173 female horseshoe crabs:
Check if the width of female’s back can explain the number of satellites attached by fitting a Poisson regression model with width.
Lets import the datasset, fit a poisson loglinear model, plot the fit and perfom a model check :
crabs = read.csv("data/crab.csv",header=TRUE)
# plot:
plot(crabs$W,crabs$Sa,col="blue4",
ylab = "Number of satellites", xlab = "width of female's back")
abline(h=c(0),col="light blue")# fit
fit.glm = glm(Sa~W,data=crabs,family=poisson)
library(gamlss)
# plot fit:
plot(crabs$W,crabs$Sa,col="blue4",
ylab = "Number of satellites", xlab = "width of female's back")
abline(h=c(0),col="light blue")
axe.x = seq(15,40,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])
lines(axe.x,f.x,col="pink2",lwd=2)## GAMLSS-RS iteration 1: Global Deviance = 923.1762
## GAMLSS-RS iteration 2: Global Deviance = 923.1762
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = -0.1836972
## variance = 2.855978
## coef. of skewness = 0.5993409
## coef. of kurtosis = 2.564461
## Filliben correlation coefficient = 0.9717529
## ******************************************************************
## GAMLSS-RS iteration 1: Global Deviance = 923.1762
## GAMLSS-RS iteration 2: Global Deviance = 923.1762
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = -0.2039895
## variance = 2.938738
## coef. of skewness = 0.5165437
## coef. of kurtosis = 2.603542
## Filliben correlation coefficient = 0.9819214
## ******************************************************************
## GAMLSS-RS iteration 1: Global Deviance = 923.1762
## GAMLSS-RS iteration 2: Global Deviance = 923.1762
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = -0.1817619
## variance = 2.766722
## coef. of skewness = 0.566886
## coef. of kurtosis = 2.728251
## Filliben correlation coefficient = 0.9772948
## ******************************************************************
# confirm lack of fit -> bin the estimates
# 2 alternative models
plot(gamlss(Sa~W,data=crabs,family=ZIP))## GAMLSS-RS iteration 1: Global Deviance = 766.424
## GAMLSS-RS iteration 2: Global Deviance = 760.0696
## GAMLSS-RS iteration 3: Global Deviance = 760.0688
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = 0.03640043
## variance = 1.195172
## coef. of skewness = 0.6257093
## coef. of kurtosis = 4.205008
## Filliben correlation coefficient = 0.9849825
## ******************************************************************
## GAMLSS-RS iteration 1: Global Deviance = 751.2934
## GAMLSS-RS iteration 2: Global Deviance = 751.291
## GAMLSS-RS iteration 3: Global Deviance = 751.291
## ******************************************************************
## Summary of the Randomised Quantile Residuals
## mean = -0.01039779
## variance = 1.043264
## coef. of skewness = -0.1133986
## coef. of kurtosis = 2.291098
## Filliben correlation coefficient = 0.9929856
## ******************************************************************
Reasonably, there is a lack of fit> the estimates are not to be trusted.